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CI,- ABSTRACT 

6 ■ 

^ ■ We investigate the influence of coherent structures on particle acceleration 

in the strongly turbulent solar corona. By randomizing the Fourier phases of a 
pseudo-spectral simulation of isotropic MHD turbulence (Re ~ 300), and tracing 

^ ■ collisionless test protons in both the exact-MHD and phase-randomized fields, it 

is found that the phase correlations enhance the acceleration efficiency during the 
first adiabatic stage of the acceleration process. The underlying physical mech- 
anism is identified as the dynamical MHD alignment of the magnetic field with 
the electric current, which favours parallel (resistive) electric fields responsible 
for initial injection. Conversely, the alignment of the magnetic field with the 
bulk velocity weakens the acceleration by convective electric fields — u x b at a 
non-adiabatic stage of the acceleration process. We point out that non-physical 
parallel electric fields in random-phase turbulence proxies lead to artificial accel- 
eration, and that the dynamical MHD alignment can be taken into account on 
the level of the joint two-point function of the magnetic and electric fields, and 
is therefore amenable to Fokker-Planck descriptions of stochastic acceleration. 

Subject headings: acceleration of particles — MHD — methods: numerical 
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1. INTRODUCTION 

Turbulent electromagnetic fields can act as particle accelerators in astrophysical situ- 
ations such as solar flares (Miller et al. 1997; Benz 2003). The turbulence is governed by 
non-linear fluid equations, and the turbulent fields therefore deviate from a simple Gaussian 
field (Adler 1981), where the Fourier phases are at random, and the field is completely spec- 
ified by its two-point function. In real space, the turbulent phase correlations manifest in 
coherent structures, with alignment of magnetic and velocity fields (Pouquet et al. 1986), 
and in non-Gaussian higher-order correlations (McComb 1990; Falkovich et al. 2001). 

The latter are not captured by traditional stochastic acceleration theory, which invokes 
the two-point functions of the magnetic field to predict the particle diffusion coefficients (e.g., 
Urch (1977); Achterberg (1981); Karimabadi et al. (1990); Urch (1991); Bykov and Toptygin 
(1993); Schlickeiser and Miller (1998); Schlickeiser (2002)). On the numerical side, some 
studies (e.g., Mace et al. 2000; Paesold et al. 2003; Arzner & Vlahos 2004) have worked with 
random-phase turbulence proxies and therefore neglected coherent structures, while others 
(e.g., Matthaeus et al. 1984; Dmitruk et al. 2003) used simulated MHD fields and took them 
into account. 

The presence of coherent structures may have conflicting effects on the acceleration 
efficiency (Dmitruk et al. 2003). On the one hand, the coherent structures may act as 
traps and detain particles on their way to high energies. On the other hand, coherent 
structures may host large electric fields which help acceleration. Which of the two effects 
prevails is hard to predict on theoretical grounds and depends on the topology and size of the 
coherent structures. A numerical exploration is thus needed to capture the whole diversity of 
possible orbit behaviour. The present article reports on relativistic test particle simulations 
in turbulent MHD electromagnetic fields, following a similar approach as Dmitruk et al. 
(2003). By applying an artificial phase randomization to the electromagnetic fields we then 
study the effect of coherent structures on test particle acceleration, and compare it with the 
Gaussian field limit. 

The astrophysical realization envisaged by the present simulation is proton acceleration 
in solar flares, with emphasis on the high-energy regime where adiabaticity of the orbit 
breaks down. Scaling to high-energy electrons or heavier ions is in principle possible, but 
not explicitly discussed. The paper is organized as follows. Section 2 describes the MHD 
fields and their randomization, and Section 3 describes the numerical generation of particle 
orbits and the physical approximations made in doing so. Section 4 presents the particle 
results, followed by a discussion of the acceleration mechanism (Sect. 5). Section 6 contains 
a summary of the simulation results, and a discussion of their validity and implications. 



- 3- 



2. ELECTROMAGNETIC FIELDS 

2.1. Random Versus Coherent Fields 

As indicated in the introduction, the focus of this paper is on the role of turbulent 
coherent structures on particle acceleration. Here we assume that particles are only subject 
to an electric field e and a magnetic field b through the (relativistic) Lorentz force (we will 
refer to a couple e and b as an "accelerator"). The electric field is taken according to MHD, 

e = -uxh + ri), Vxb = j, (1) 

where u is the bulk velocity, j is the electric current and rj is the magnetic diffusivity. In 
Eq. (1) and from now on we use Alfven units, i.e., the magnetic field is scaled to the 
Alfven velocity, b = B/^//z p, and measured in the same units as the bulk velocity. Another 
assumption of the present study is that the particles' motion take place in a cubic domain 
with periodic boundary conditions. This implies that at any instant the MHD fields can be 
viewed as superpositions of waves and decomposed in Fourier modes: 

u(x) = $>(k)e lk - x , (2) 

b(x) = ^b(k)e lkoc . (3) 

In the above expansions, u(k) and b(k) are complex numbers which satisfy, u(— k) = u*(k), 
b(— k) = b (k) and can be written as, 

u(k) = u r (k)e ia(k) , b(k) = b r (k)e i/3(k) (4) 

with real-valued u r , b r , a and (3. Constructing a numerical field like (4) thus requires 
information both on the modulus and the phases of its Fourier modes. In homogeneous 
and isotropic turbulence, the modulus of Fourier modes are usually set to given values 
that depend only on the norm of the corresponding wave vector to mimic a given energy 
spectra. Assigning appropriate phases to the Fourier modes is a much more complex task. 
In a real turbulent field, these phases have to be correlated in such a way as to produce a 
field with subtle and intricate structures. Prescribing them numerically represents the same 
complexity as prescribing the electromagnetic field at each point in space and is thus not 
practical. The most simple way of circumventing this problem is to assign random numbers 
in the interval [—it, 7r] to the phases and impose a given distribution, most often the uniform 
distribution to respect isotropy. The fields generated can then be qualified as stochastic but 
still have a prescribed energy spectra. More realistic phases can be obtained by generating 
the electromagnetic field through direct numerical simulations (DNS) of the MHD equations, 

d t Ui = -diip/p) - ujdjUi + bjdjbi + uAui, (5) 
dA = -ujdjbi + bjdjUi + rjAbi, (6) 
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where p is the sum of the kinematic and magnetic pressures and v is the kinematic viscosity. 
The magnetic field contains no mean contribution. Since we assume incompressibility, p is 
obtained by solving a Poisson equation and is not an independent variable (e.g., McComb 
1990). Although the initial condition used in (5) and (6) consists of fields with random 
phases, it is observed that after some time, the resulting electromagnetic field exhibits co- 
herency that deviates strongly from its random phases counterpart (e.g., Biskamp & Miiller 
1999). Although the DNS approach is very appealing, it should be stressed that it is very 
costly and that the values of Reynolds numbers observed in typical turbulent systems are 
out of reach of present-day supercomputers. In that sense, the fields obtained through DNS 
also constitute models for real-life turbulence. 

Our objective is here to study two kinds of accelerators: a random phases one and a 
coherent one obtained through DNS. To make things consistent, we construct both accel- 
erators in such a way that they have exactly the same spectra and only differ through the 
value of the phases of their Fourier modes. As a shorthand, the random phases accelerator 
will be denoted (SB)sto while the DNS accelerator will be denoted (SB)dns- in the next 
section we detail the construction of (£B)sto and (£B)dns- 

2.2. MHD Simulations 

In order to obtain a suitable DNS electromagnetic field, we solve equations (5) and 
(6) using a pseudo-spectral code in a cubic domain of side length l x — l y — l z — 2ix and 
impose periodic boundary conditions. Our simulation code is dealiased and uses a 4th-order 
accurate Runge-Kutta time advancement scheme. As explained above, the initial condition 
for u,i and hi are random phases fields of the form (4) whose spectra match the one observed in 
the Comte-Bellot-Corrsin grid turbulence experiment at stage 2 (Comte-Bellot and Corrsin 
1971) (this choice of spectra is motivated by the fact that the corresponding turbulence can 
be resolved with a resolution of 256 3 Fourier modes). The only extra requirement imposed 
on the Fourier modes is incompressibility, k-b(k) = k-u(k) = 0. As time advances, the flow 
becomes "physical" and we stop the simulation when the skewness of the velocity derivative 
dui/dxi reaches a quasi constant value (Batchelor 1982; McComb 1990; Gotoh et al. 2002; 
Verma 2004). This criteria is standard and indicates that turbulence has had time to develop. 
Here, the initialisation phase has taken about 2.5^ where is the Alfven time defined by 
tA = L/va, with L the (magnetic or velocity) integral scale-length and v\ = (|b| 2 )/3 (from 
here on, angular brackets (...) denote a spatial average). Note that va is determined by 
the (rms) magnetic fluctuations and not by a background magnetic field. Some of the 
characteristics of the flow at 2Ma are listed in Table 1. There is approximate equipartition 
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between |b| 2 and |u| 2 , and both dissipate at the same rate (magnetic Prandtl number ~ 1), 
so that the magnetic and velocity characteristic scales are similar. The Reynolds number and 
magnetic Reynolds number are about 300. The magnetic field extracted from the DNS at 
this final time and the electric field obtained through (1) constitute the accelerator (SB)dns- 
To construct (£B)sto we modify the phases of b(k) and u(k) according to, 



Resolution 256 3 

Box size (l x x l y x l z ) 2tt x 2n x 2n 

Rms velocity u = yj 2.20 

Rms magnetic field b = \J 2.39 

Viscosity v 0.006 

Diffusivity rj 0.006 

Kinetic dissipation e u 15.28 

Magnetic dissipation e h 20.57 

Integral length-scale of u (L u = 37r/4 x (J k~ 1 E u (k)cIk/ J E u (k)cIk)) 0.79 

Integral length-scale of b (L& = 37r/4 x (J K^ 1 Eb(n)d,K/ f Eb(K)dK)) 0.71 

Longitudinal Taylor microscale of u, A 2 = 15uu 2 /e u 0.169 

Longitudinal Taylor microscale of b, A 2 = lhrjb 2 /e h 0.158 

Eddy turnover time (E u /e u ) 0.421 



Table 1: Turbulence characteristics of the velocity and magnetic fields used to construct 
(SB)dns- E u (k) and E^k) designate the velocity and magnetic field spectra, e.g., E u = 
|m 2 = f E u (k)cIk. Angular brackets denote a spatial average. 



a l (k)^«,(k) + 0(k), A(k)^A(k)+v?(k), (7) 

where 0(k) and y?(k) are random numbers distributed uniformly in the interval [— n, ir] and 
which satisfy 0(k) = — 0(— k) and <p(k) = — ip(— k). The electric field is then computed from 
Eq. (1). The above transformation ensures that 1) the phases of the new fields (SB)sto 
contain a highly random part; 2) the spectra of u and b remain unchanged; 3) the divergence 
of b remains equal to zero; 4) the ideal MHD condition e ■ b = is preserved if 77 = 0, so 
that the randomization does not introduce parallel electric fields. 



2.3. Intrinsic Properties of the Accelerators (£B)dns and (£B)sto 



The whole purpose of building an accelerator using DNS of the MHD equations is to 
obtain an electromagnetic configuration that has more characteristic features of turbulence 
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than only the correct energy spectra. In this section we describe some of the differences 
between (£B)dns an d (£B)sto that are relevant from the point of view of turbulence and 
particle acceleration. See Verma (2004) for a recent review on MHD turbulence. 

Figure 1 presents 3D views of the magnetic (|b| 2 ) and electric (|e| 2 ) energy densities. 
Although this corresponds only to a visual impression, it appears clear from the graphs that 
the electromagnetic field composing (£B)ons hosts more coherent, filamentary structures 
than the one composing (£B)sto- 111 the latter, structures appear noisier and span regions 
of smaller extent. 
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Fig. 1. — Magnetic (top line) and electric (bottom line) energy densities of the (£B)dns 
(left column) and (£B)sto (right column) fields, similar as in Dmitruk et al. (2003). Blue 
and red regions indicate respectively low and high values. 
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Also, it appears from Fig. 1 that the electric field has more extreme values (dark 
red) than the magnetic field, as already noticed by Dmitruk et al. (2003). More quantitative 
information on this can be obtained from the distributions (one-point functions) of individual 
field components. In the case of (SB)sto it is observed, as expected, that U\ and b± have 
a nearly Gaussian distribution (other components behave similarly). The distribution of e± 
(again for (£B)sto) is very close to a zero-th order modified Bessel function of the second 
kind, which decays exponentially for large values of its argument. This is again expected, 
since at high Reynolds number e x is dominated by the term — u x b which behaves as the 
product of two Gaussian variates for (£B)sto- In the case of (£B)dns, the distribution of 
the field components are similar to their random-phase equivalents. Only for large values do 
they exhibit moderately pronounced tails resulting from higher intermittency. 
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Fig. 2. — Histograms for the longitudinal derivatives of magnetic (left) and electric (right) 
field components. Solid - DNS; dashed - random phase. 
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A better measure of coherency is contained in the distributions of the longitudinal 
derivatives such as dibi and die x (cf. Gotoh et al. 2002) (Fig. 2). Again, the 2nd and 
3rd components behave similarly and are therefore not shown. In the DNS case, both d\bi 
and d\e\ depart significantly from their random-phase equivalents. The distributions exhibit 
a large number of extreme values (tails in the distribution). These events allow coherent 
structures to be present in the signal since they correspond to regions of space where sharp 
transitions "shaping" the signal can occur. The strong differences in this diagnostic between 
the DNS and random phases case reflect the marked visual differences observed in Fig. I. 
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Fig. 3. — Diamonds: the excess of the kurtosis of A6; = 6j(x + rej) — 6j(x) between DNS and 
STO fields, averaged over % = 1,2,3. Inlets: some representative histograms of Abi (average 
over i=l,2,3); black and gray lines represent DNS and STO cases. See text. 




Fig. 4. — Distributions of the angle cosines between b and u (top), and b and j (bottom). 
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The coherent structures also manifest in the longitudinal structure functions (Kol- 
mogoroff 1941a,b; McComb 1990; Frisch 1995; Sorriso-Valvo et al. 2000; Falkovich et al. 
2001; Gotoh et al. 2002; Cho et al. 2002) which are a natural extension of the longitudinal 
derivatives discussed above. The magnetic longitudinal structure functions are defined by 
D p (r) = (|[b(r + x) — b(x)] • f| p ) and measure the variation of the p-th moment of the 
longitudinal magnetic field increment with distance r (hat denotes a unit vector). Figure 3 
(inlet) shows the distribution of these increments along fj, averaged over i = 1,2,3. With 
increasing r, the increments become decorrelated and their distribution becomes Gaussian. 
The transition can be monitored by means of the kurtosis n(r) = D A (r) / D 2 {r) 2 \ Fig. 3 (dia- 
monds) shows the DNS-to-STO kurtosis ratio as a function of r. There is a clear systematic 
excess of the DNS kurtosis at small r, and one may infer from Fig. 3 that on scales r > 1, 
the DNS and STO fields behave similarly. In particular, the DNS field becomes independent 
at separations which are long compared to the integral length-scale of u and b (Tab. 1). 

Another distinction between (£B)dns and (£B)sto, which will prove helpful in under- 
standing the particle acceleration, is the relative orientation of b with respect to u and j. 
The nonlinear MHD equations (6) and (5) tend to align b with j in order to minimize the 
j x b force. The (approximate, Ret, = 283) frozen-in condition tends, then, to align b with 
u (Biskamp 1993; Dar 1998; Debliquy et al. 2005). This preference for alignment is demon- 
strated in Fig. 4, showing the distribution of the orientation cosines b • u and b • j within the 
simulation cube. The DNS distributions (black line) accumulate at the extreme values ±1, 
corresponding to ±180°, more pronounced for b • j than for b ■ u. After the phase random- 
ization is applied (gray line), this accumulation vanishes and the distribution of b • j (botom 
panel) become isotropic. The residual anti-correlation between b and u (top panel, gray 
line), corresponding to a negative cross-helicity, is dominated by few small- |k| modes, and is 
within statistical errors. Note that for the mixed-field correlations, the coherent structures 
already manifest at the level of b-u and bj, whereas one has to go to higher order moments 
if b and u are separately considered (Fig. 3). 



3. PARTICLE ORBITS 

3.1. Physical Environment 

An ensemble of collisionless charged particles is amenable to the MHD description if 
their orbits can be described by the drift velocity e x b/|b| 2 + vy, where '||' refers to the 
local magnetic field. The drift velocity then equals the MHD bulk velocity u. This de- 
scription requires that the Larmor radius is much smaller than the radius of curvature of 
the magnetic field, and that v\\ is smaller than the Alfven velocity (e.g., Littlejohn (1982); 
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Biichner & Zelenyi (1998)). In the quiet solar corona, the ordering v < va holds true for 
both electrons and ions. During solar flares, when turbulence arises, a small fraction of the 
particles reaches relativistic energies, and it is this high-energy population which is treated 
here as test particles. The non-energetic bulk provides the MHD acceleration environment. 
The acceleration to relativistic energies is known to be rather rapid. In impulsive flares, hard 
X ray signatures of relativistic electrons appear within fractions of a second, and Gamma 
ray signatures of relativistic protons - if any - appear within a few seconds (Miller et al. 
1997). The MHD structures, in contrast, evolve on a much longer time scale. 

The solar corona, in which acceleration is thought to occur, is a relatively collisionless 
environment. Assuming a magnetic field strength of 10~ 2 T, the electron and proton gyro 
times are in the order of nanoseconds and microseconds, respectively. In contrast, the ion- 
electron collision time is in the order of tens of milliseconds if we adopt a typical temperature 
T e = Ti = 10 6 K and density 10 16 m~ 3 . The particles of the background plasma thus perform 
many gyrations before a Coulomb collision is encountered, and this becomes even more 
true for high-energy test particles, because the collision rate rapidly decreases with energy 
(Dreicer 1959; Helander & Sigmar 2002). While the background plasma relaxes towards a 
Maxwellian on time scales of tens of milliseconds, a relativistic test proton remains virtually 
collisionless. 

For particles which move much faster than the background plasma, the MHD fields may 
be regarded as static. This requires v ^> va- In the solar corona, va is in the order of 1000 
km/s, and the condition v ^> v a is fulfilled for relativistic particles. However, acceleration 
begins at Alfvenic velocities, and the question rises whether a static approximation applies 
from the very beginning. This depends on the time spent at low velocities. As outlined 
above, this time is observationally known to be short compared to the MHD time scale in 
the real solar corona, and in our orbit simulations we will find that particles reach rela- 
tivistic velocities within some 0.01^- The accelerated particles therefore spend most of the 
simulation in a state with v >> va, and we set 



throughout the whole orbit simulation. The stationarity assumption (8) admits conservation 
laws which facilitate the numerical benchmarking and the physical understanding of the 
acceleration process. In this sense, Eq. (8) also represents a simplifying working hypothesis. 

The test particles orbits are measured in the same units as the MHD fields. The dimen- 
sionless equations of motion read 





~dt 



= v 



(9) 
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-L- = a(e + vxb), (10) 

where 7 = (1 — v 2 /c 2 )~ 1 / 2 is the Lorentz factor, and the effective charge a = VLt A defines the 
particle time scales relative to the MHD time scale (Dmitruk et al. 2003), with VL the proper 
rms gyro frequency. Assuming a coronal magnetic field of 10 _2 T and a density of 10 16 m~ 3 , 
we have va = 2182 km/s = 0.0073c (thus c=137.61 in numerical units), and focusing on 
protons as test particles, we set a = 10 3 . For display purposes, time is measured in units 
of the (dimensionless) proper gyro time tc = 27t/(o;a/ (|b| 2 ) = 3.6/ (av a), corresponding to 
about 5 ■ 10~ 6 s in real time. All test particles start out at random position and in random 
direction with initial velocity v = 0.2v A , representative of protons with T = 10 7 K. This 
ensures that the test protons can be considered as approximately collisionless, and that the 
initial Larmor radius is smaller than the numerical cell. The orbits are exactly integrated, 
without recourse to gyrokinetic approximations. This approach benefits from simplicity and 
rigor, but allows only relatively small populations to be simulated over relatively short times. 
See Section 6 for further discussion of the limitations of the method. 



3.2. Numerical Implementation and Benchmarking 

Equations (9) and (10) are integrated by a 4th order Runge-Kutta scheme with constant 
time step. The electromagnetic fields are taken from a snapshot of the pseudo-spectral MHD 
simulation, transformed to real space, and linearly interpolated to the actual particle posi- 
tions. A parallel (distributed memory) implementation is used to enhance the computational 
power. The diagnostics includes energy histograms and the second moments of momentum 
and space displacements. 
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Fig. 5. — Benchmarking of energy conservation (Eq. 11). Top: relative error along an 
example orbit. Bottom: ensemble average, with the shaded region representing the ±5cx 
deviations. Time is measured in units of the proper rms gyro time tc- 
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The static assumption (Eq. 8) implies that the vector potential evolves linearly with 
time, and that the scalar potential is independent of time. Therefore the identity 



must hold, where e so i is the solenoidal (V • e so i = 0) part of the electric field and H Q is 
the initial energy. Numerically, and e sol are computed from the Fourier amplitudes e(k) 
by decomposition into parallel and perpendicular components with respect to k. To make 
the decomposition unique, the potential is gauged such that (0) = 0. Equation (11) can 
be used to assess the numerical accuracy and to choose an appropriate time step. To this 
end, the right hand side of Equation (11) is numerically integrated along the trajectories 
and compared to the lefthand side (Fig. 5). The observed discrepancy AE depends on the 
time step At and on the smoothness of the MHD fields; both small At and smooth fields 
help to fulfill Equation (11). In view of the fact that real acceleration increases the particles' 
kinetic energy by orders of magnitudes, an rms discrepancy of 1% on the 2a level between 
the left- and righthand sides of Equation (11) at the end of the simulation is considered as 
acceptable. This is achieved if At < 0.01 x (cnu) -1 - 

The Runge-Kutta integration scheme was also tested against a simpler leapfrog scheme 
and it was shown that it provided a significant increase in precision. 




(11) 
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4. RESULTS 

In this Section we summarize the outcome of the particle simulations, focusing on the 
second moments and correlation functions of the dynamical quantities (x, p). The average 
over particles will be denoted by ((...)) in order to distinguish it from the spatial average (...). 
The particle ensembles comprise 50.000 particles. 
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Fig. 7. — Mean square displacement in correlated (black) and decorrelated (gray) fields. The 
dashed and dotted lines indicate linear and quadratic growth. 



4.1. Real-Space Diffusion 



A first qualitative difference between the (£B)dns and (SB)sto cases concerns the 
intermittency of their real-space trajectories. Figure 6 shows sample orbits in (SB)dns 
(left) and (£B)sto (right) fields for times up to t/tc < 1.5 • 10 3 . As can be noticed, the 
(SB) dns orbits tend to have long sojourns between clusters where particles are temporarily 
trapped, and depart more from their origin than the (£B)sto trajectories. This tendency, 
observed for the few samples in Figure 6, reflects in a full population average of the mean 
square displacement (Fig. 7). Here we see that the (£B)dns average grows somewhat faster 
than the (£B)sto average during the first 200 gyro times, from where on a relative excess 
of about 30% is preserved up to t/tc ~ 10 4 . During this time interval, the mean-square 
displacement grows approximately linearly with time (dashed). Later on, the (£B)sto curve 
overtakes the (£B)dns curve, and grows again faster than linearly with time. The bends 
observed in Fig. 7 indicate changes of acceleration regimes which will be discussed below. 

The distance traveled by the test particles ultimately exceeds the side length of the 
MHD simulation cube, and the particles cross many periodic continuations of this cube. 
However, the particle orbits themselves are not periodic, because opposite points like (x, y, 0) 
and (x, y, l z ) are not magnetically connected (this is a consequence of the fully 3-dimensional 
nature of the MHD simulation), and because the particles may leave their magnetic field lines. 
Rather, the particle orbits have mixing and ergodic properties. Therefore, the periodicity of 
the MHD fields is not expected to affect the particle results. 
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4.2. Momentum Diffusion 



Next we consider the diffusion in momentum space (momentum p = 7V is a better 
indicator for diffusion than the velocity because it is not limited by the speed of light). The 
evolution of the mean square momentum is shown in Figure 8 for both (£B)dns an d (£B)sto 
ensembles. The (£B) DNS curve grows faster at 10 < t/t c < 10 2 , which is the typical particle 
crossing time of the coherent structures. Later, when many coherent structures have been 
visited (t/t c ~ 5- 10 3 ), the (SB)dns curve overtakes. The average terminal (t/tc = 1-9- 10 4 ) 
energy is ((7)) = 6.4 for the DNS case and ((7)) = 8.0 for the STO case. Thus our protons 
are accelerated to relativistic energies. 
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Fig. 10.— Energy distributions at t/t c = 643, 1264, 2299, 4369, 6439, 10579, 16789 (left to 
right). 
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4.3. Momentum-Displacement Correlation 



The propagation of diffusion from momentum space (where (SB)dns and (SB)sto ac t) 
into real-space is governed by the mixed correlation -j|((|x(£)— x(0)| 2 )) = 2(([x(£) — x(0)]-v(t))). 
This quantity is always positive due to the statistical preference for particles having traveled 
away from the origin to have a velocity component in the same direction. Since the coherent 
structures are able to channel the particles along filaments, they impose a correlation between 
position and momentum, and we may expect (([x(£) — x(0)] -p(t))) to grow faster in (£B) DNS 
than in (SB)sto during the first structure-travel time. This is indeed the shown 
in Fig. 9. The plateau at intermediate times 100 < t/t c < 2000 indicates the regime of 
ordinary diffusion. 
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Fig. 11. — Particle anisotropy in the (£B)dns an d (£B)sto fields. 
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Fig. 12. — Orbital radius of curvature versus time. Dots: sample orbits; solid line: average. 
For comparison, the numerical cell 8 and the Taylor microscale A& are also indicated. 
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4.4. Energy Spectrum 

The evolution of the full energy distribution is shown in Fig. 10. The particles start 
out at t=0 with kinetic energy 0.0037mc 2 , and spread then to higher values. Note that for 
t/tc < 3 • 10 3 , the DNS fields are more efficient accelerators than the STO fields, while 
the converse holds true at later times. The spectral shape has slightly positive (negative) 
curvature for t/t c < 5 ■ 10 3 (t/t c > 5 • 10 3 ), indicating that the distribution functions decay 
sub-exponentially (super-exponentially). When fitted by a power law N(E) ~ E~^ in the 
sub-exponential range 400 < t/tc < 2000, power law indices £ = 4. 5. ..3. 7 are found. Similar 
exponents have been observed in medium-size solar flares, although for electrons, and smaller 
energies E/m e c 2 < 0.2 (Grigis and Benz 2004). 

5. ACCELERATION MECHANISM 

The ultimate cause of particle acceleration is the electric field which is now examined 
in some detail, together with the microscopic nature of the particle orbits. From a particle 
dynamics point of view, the main classification of e is into potential-solenoidal and parallel- 
perpendicular (to b) components. From an MHD point of view, one may further distinguish 
between convective and dissipative fields. 

5.1. Solenoidal and Potential Electric Fields 

By virtue of the 'snapshot' condition (Eq. 8), the potential part — V0 of the electric 
field does not contribute to acceleration on time scales which are long compared to the 
crossing time of local potential minima. However, the initially spatially uniform distribution 
of particles may become non-uniform at later times, resulting in a shift of the average kinetic 
energy which would manifest in Figs. 8 and 10. This shift is bounded by the magnitude of 
the potential fluctuations A7 = aA(f)/c 2 < 0.08, which is small compared to the energies 
reached in the tail of the energy distribution (Fig. 10). It is thus the solenoidal electric field 
e so i which is responsible for acceleration. From an MHD point of view, e so i is made up of 
the dissipative field rjj plus parts of the convective field — u x b. Numerically, the individual 
(rms) contributions are e£J c s onv = 8.0, e s s ^° onv = 9.35, and e s D ^ iss = e ^° diss = 0.35. Thus the 
solenoidal electric field is somewhat weaker in the DNS than in the STO case, in contrast to 
the enhanced DNS acceleration efficiency at t/t c < 5 ■ 10 3 . 
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5.2. Parallel and Perpendicular Electric Fields 

The picture, though, changes if we distinguish between parallel and perpendicular com- 
ponents where e™| = 0.24 > ej^jj = 0.20. Initially (t — 0), the Larmor radius is smaller 
than the size of the magnetic inhomogeneities, so that the particles are in a gyrokinetic 
regime. The ratio of the (rms) forces on the particles is |v x b| 4- |e_|_| 4 |e||| ~ 1.74- 12 4-0.3, 
so that the particles, typically, exhibit a strong e x b drift. The drift velocity v d = exb/|b| 2 
thus approaches va, and the maximal drift Larmor radius r* L = ^a\e±\(tc/'2) 2 exceeds its 
initial (thermal) value v±/ (ab) by up to an order of magnitude. This is a consequence of the 
strong turbulence assumption, where Va is given by the perturbations and not by a (dom- 
inant) background magnetic fields. Although the particles gain substantial perpendicular 
momentum during the first half gyration, this is approximately restored in the second half, 
and so forth, so that there is little net gain from ej_. The parallel field ey, however, yields 
direct parallel acceleration. As a consequence, the velocity distribution becomes anisotropic 
(Fig. 11), with a maximum anisotropy reached at t/tc ~ 50. At later times, the velocities 
isotropize. One reason for this is the (reversible) conversion of parallel into perpendicular 
momentum at converging magnetic field lines. This mechanism is supported by the observa- 
tion that ((|b| 2 )) is slightly (0.04 2 ) enhanced during the isotropization phase 30 < t/t c < 300, 
as would be expected when particles penetrate high-|b| domains. 

A second cause of the isotropization observed in Fig. 11 is the fact that the gyrokinetic 
approximation breaks down with increasing energy. In order to monitor the transition from 
gyrokinetic to eventually free orbits we compare the orbital radius of curvature p = 1/K 
to the characteristic scales of the magnetic field, where K 2 = (|x| 2 |x| 2 — |x ■ x| 2 )/|x| 3 . The 
evolution of p is shown in Fig. 12 for a sample of 50 particles (dots), with the average {(p)) 
marked by solid line. Initially, ((p)) is smaller than the Taylor microscale A& and also smaller 
than the numerical cell S (the electromagnetic fields are interpolated to the actual particle 
position). It then rapidly grows during < t/tc < 50, where it approaches the cell size and 
the growth slows down. The quantity r* L defined in the foregoing paragraph is found equal 
to the maximum of p during the first gyration within ±10%. Based on Fig. 12 and the 
criterion ((p)) > 5 one may infer that, on average, the gyrokinetic regime is left after a few 
hundred gyro times. At this point the character of the orbit and the acceleration mechanism 
change. 

A picture thus arises where acceleration is initiated by the parallel component of the 
solenoidal electric field ('injection'), and is then taken over by the full non-potential electric 
field. While the parallel component is somewhat larger in (SB)dns, the total solenoidal 
field is larger in (£B)sto- Thus the DNS fields are initially stronger accelerators, but the 
STO fields become more efficient at a later (non-adiabatic) stage. In both cases, the parallel 
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electric field is entirely dissipative. The initial (parallel) acceleration phase corresponds to 
the first phase of Fig. 7 (t/t c < 200). 

The fact that the parallel field is weakened by the phase randomization, while the 
perpendicular electric field is strengthened, is easily understood from Fig. 4. The non- 
linear MHD alignment of b with j favours small this preference is lost in the phase 
randomization so that increases. Similarly, the MHD alignment of u and b is lost 
in the phase randomization, and |u x b| increases. Since the convective field numerically 
dominates the dissipative one, the solenoidal electric field as a whole increases when the 
phase randomization is applied to u and b. 
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Fig. 13. — Similar to Fig. 8, but with e directly phase randomized (instead of computing it 
from the phase randomized b and u). The time axis spans a smaller range than in Fig. 8. 
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5.3. Direct Phase Randomization of E 

The crucial role of parallel electric fields can also be demonstrated by an alternative 
simulation setup, where e is directly phase-randomized rather than computed from the phase- 
randomized u and b. The ideal MHD condition e • b = is then broken even in the absence 
of resistivity, while the power spectral densities of b and e remain unchanged. We may 
thus ask what happens to the particles in such a directly phase-randomized electric field. 
This procedure has a certain interest on its own because of its implication on random-phase 
Fokker-Planck descriptions of stochastic acceleration. 

It turns out that the direct phase randomization of e leads to a dramatic enhancement 
of acceleration. Figure 13 shows the evolution of (( |p(t) — p(0)| 2 )) for the DNS and directly 
phase-randomized cases. The curve corresponding to direct phase randomization of e is 
labeled STO (E). As can be seen, the DNS and STO (E) curves separate after about half 
a gyro period, and the direct phase randomization of the electric field yields energies which 
are two orders of magnitude above the (SB)sto results. This dramatic enhancement of 
acceleration is entirely due to (unphysical) parallel electric fields. The local maximum of 
the DNS curve at tc*/2 marks the first half gyro phase, where the velocities are, on average, 
reversed. At this point, the particles in the directly phase-randomized electric field have 
already gained super- Alfvenic velocities so that the effect of gyration is no longer visible. 

6. SUMMARY AND DISCUSSION 

We have performed a numerical study of test particle acceleration in strong MHD turbu- 
lence (b, u) and its phase-randomized version, with the electric field given bye = — uxb+r^j. 
It is found that the dynamical alignment of b with j and u in true MHD affects the division of 
the electric field into parallel and perpendicular components, in the sense that parallel (resis- 
tive) components are favoured and perpendicular (convective) components are disfavoured. 
As a consequence, the true MHD turbulence provides more efficient direct acceleration of 
super-Dreicer particles than its phase-randomized version. If the particles reach sufficiently 
high energies to leave the gyrokinetic regime, so that they can become accelerated by the 
solenoidal part of the u x b force, the dynamical alignment of u with b yields weaker accel- 
eration than predicted from the phase-randomized fields. For collisionless protons in MHD 
turbulence with magnetic Reynolds number ~ 300, this is the case for protons after a few 
1000 gyro times. Numerically, the acceleration efficiencies of the true MHD and phase- 
randomized fields differ by factors of order two. It is expected that the discrepancy would be 
larger for super-Dreicer electrons, since these remain longer in the adiabatic regime where 
the enhancing effect of (b,j) alignment is not yet competed with the diminishing effect of 
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(b, u) alignment. 

With regard to solar flares envisaged as a astrophysical application, a few words of 
caution are appropriate. The present simulation involves only dimensionless quantities such 
as tl/Xij or tl/5. In order to fix the absolute physical scaling, and since we were interested 
in the possibility of non-adiabatic acceleration, we have chosen protons as test particle. 
Using typical values of the coronal magnetic field strength, and taking the resolution of 
the simulation into account, the protons become non-adiabatic (tl > 8) at v > 0.6 c, and 
reach energies as high as 10 GeV (10 proton rest masses). Such high energies are commonly 
not observed; the highest energies from in-situ observations in the interplanetary space are 
about 100 MeV/nucleon (Reames et al. 1992), and proton energies inferred from pion decay 
photons reach up to a few 100 MeV/nucleon and occasionally up to a few GeV (Kanbach et al. 
1993). (Only in the very largest flares such as of October 28, 2003, ground-based detections 
of relativistic protons and neutrons show energies up to 10 GeV (Bieber et al. 2004).) Also, 
the simulated acceleration is very rapid; assuming (|b| 2 ) ~ (100G) 2 , the simulation duration 
is only 0.1s in real time. There are two main reasons for the violent acceleration in the 
present simulation: the neglect of collisional and radiative losses, and the freedom in the 
absolute scaling of the simulation system. The latter is chosen such that one numerical cell 
is about 200 initial (|v| = v ) Larmor radii; this corresponds to 5 ~ 100m in real space. The 
resulting steep gradients - which are not detectable by present-day observations -, together 
with a larger-than-coronal resistivity (i?e& ~ 300), yield highly efficient accelerators, which 
allow us to rigorously follow the particles to relativistic energies before becoming swamped by 
accumulated numerical errors (and running out of computing time). The simulation system 
is thus affected by computational considerations, and presumably over-estimates the absolute 
coronal acceleration rate. However, the relative influence of dynamical MHD alignment on 
particle acceleration should be correctly reproduced. 

Another caveat stems from the test particle nature of the present simulation, which 
does not allow to address the question of absolute numbers of accelerated particles. The 
test particles are considered as resulting from a runaway process which is out of the scope 
of this simulation, so that we cannot directly relate their number to the number of thermal 
background particles. However, it is clear that the test particles must not gain more energy 
than is available in the MHD fields. This sets an upper limit on their number, and also 
on the duration of the simulation. (As the simulation does not include loss processes, the 
test particles can gain arbitrary large energy.) In our simulations, the test protons pick up 
Alfvenic velocities within one gyration. This brings them to kinetic energies of the ambient 
(low-beta) coronal plasma, which is assumed to be in equipartition, (|u| 2 ) = (|b| 2 ). In the 
course of the simulation, the kinetic energy of the test particles increases by a factor 10 4 . 
This should not exceed a small fraction (say, e = 10~ 2 ) of the magnetic field energy, such 
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as not to perturb the MHD configuration. Therefore the test particles represent a fraction 
/ = 10~ 6 of all particles only. This fraction is small, but not unreasonably small for runaway 
protons, which must be substantially faster than the electrons in order for runaway to occur 
(Dreicer 1960). It is also in agreement withe the observation that most solar flares do not 
show signature of high-energy protons (but only of high-energy electrons). However, if the 
test particles are to represent more than 10 -6 of all particles, then the acceleration should 
be stopped at earlier times. Since it is found that the test particle energy increases roughly 
linearly with time, the maximum simulation duration is t max /t c ~ e/f. In fact, it might 
happen that the MHD description brakes down after a few 100 gyro times, so that the 
non-adiabatic regime where the (SB)sto fields become more efficient accelerators is never 
reached. 

Returning to the more technical aspects, and motivated by the findings of Sect. 5.3, we 
would like to stress that violation of the ideal MHD condition b • e = by random-phase 
approximations unavoidably leads to artificial particle acceleration. If, in theory or numer- 
ics, a turbulence proxy is generated by superimposing Alfven modes (5h^, £e k )e* k ' x_;i ( k ' bo * ) * 
with <5efc = b x 5bk, then a parallel electric field occurs which is of second order in 8b. 
For example, if the 5h^ are random isotropic in the plane perpendicular to b , with ori- 
entation uncorrelated to k, then ((e • b) 2 ) = jbl((5b) 2 ) 2 ; for 5b <C b , the parallel electric 
field has thus standard deviation ^((5b) 2 }. Although this is small compared to the per- 
pendicular electric field (e\\/e± ~ 5b/b ), it may give rise to secular growth of the parallel 
particle momentum. Therefore, stochastic acceleration in turbulence proxies constructed 
from superimposed Alfven waves might in fact over-estimate the acceleration in true MHD 
turbulence. 

Finally, we point out that, although the coherent structures do not manifest in the two- 
point functions (or power spectral densities) of b and e alone, they contribute to their mixed 
correlations. This fact can, in principle, be taken into account in Langevin- or Fokker-Planck 
formulations of stochastic acceleration (e.g., Brissaud & Frisch 1974; Van Kampen 1981; 
Risken 1989; Urch (1977, 1991); Gardiner 1997; Schlickeiser and Miller (1998); Schlickeiser 
(2002); Falkovich et al. 2001), which deal with the transfer of spatial correlators of the elec- 
tromagnetic fields into the particle diffusion coefficients. Given the correlations (&i(0) bj(r)), 
(bi(0) ej(r)), and (ej(0) e, (r)), one may compute, for example, the momentum diffusion coeffi- 
cient D pp = J °°((5f(0) -5f(r))) dr along the drifting gyrocentre x (r) = x(0) + fj|(0)br + v d r 
(Sect. 5.2), where 5f = e(x) + v x b(x) for strong turbulence. The quantity D pp may 
then be evaluated numerically, as follows. First we choose a large (10 7 ) number of random 
initial data (x(0),v(0)) with uniform position x(0) and isotropic velocity of fixed modulus 
|v(0)| = v o = 0.2va, like for the simulated test particles. Then, for each (x(0),v(0)) the 
instantaneous (straight) gyrocentre orbit x (r) is computed, and the electromagnetic force 
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Sf(r) is probed along this orbit. The average of ({Sf(0) ■ Sf(r))) over all particles (initial 
data) is performed, and the integral in the definition of D pp is approximated by a sum of 
sufficiently resolved increments. Proceeding this way we obtain D pp = 0.03 for (SB)ons 
and D pp = 0.018 for (SB)sto- The ratio of these values reflects the excess of ((p 2 ))dns over 
(0° 2 ))sto observed in Fig. 8 for t/tc < 1000. A detailed study of the modification of D pp by 
the cross-correlation (frj(0) e_,-(r)) has, to the authors knowledge, not yet been carried out. 

The authors thank R. Kallenbach and M. Verma for helpful discussions. L.V. is sup- 
ported in part by the Research Training Network (RTN) 'Theory, Observation and Sim- 
ulation of Turbulence in Space Plasmas', funded by the European Commission (contract 
No.HPRN-eT-2001-00310). 
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